/* ****************************************************** */
/* * Borwein Method for finding PI - Sequential Version * */
/* ****************************************************** */

#include <gmp.h>
#include <stdio.h>
#include <stdlib.h>

/* Number of iterations */
#define NUM_ITERS 2000
/* Size of array. For best results, NUM_ITERS % ARRAY_SIZE should be 0 */
#define ARRAY_SIZE 2000

/* Global array of mpf_t variables */
mpf_t a[ARRAY_SIZE+1], aSqrt[ARRAY_SIZE+1], b[ARRAY_SIZE+1], pn, pnp;

int main(int argc, char **argv) {
	int k, i;
	mpf_t bBottom, pBottom;

	/* Define precision */
	if(2 == argc)
		mpf_set_default_prec(atoi(argv[1]));
	else if(argc > 2) {
		printf("Usage: %s [prec]\n", argv[0]);
		exit(1);
	}

	/* Set a[0] = sqrt(2) */
	mpf_init_set_d(a[0], 2.0);
	mpf_sqrt(a[0], a[0]);
	/* Set sqrt(a[0]) */
	mpf_init_set(aSqrt[0], a[0]);
	mpf_sqrt(aSqrt[0], aSqrt[0]);
	/* Set b[0] = 0 */
	mpf_init(b[0]);
	/* Set p[i] = p[0] = 2 + sqrt(2) */
	mpf_init_set_d(pn, 2.0);
	mpf_add(pn, pn, a[0]);
	/* Initialise p[i+1] */
	mpf_init(pnp);
	/* Initialise auxiliary variables */
	mpf_init(bBottom);
	mpf_init(pBottom);

	/* Initialise the rest */
	for(i = 1; i <= ARRAY_SIZE; i++) {
		mpf_init(a[i]);
		mpf_init(aSqrt[i]);
		mpf_init(b[i]);
	}

	for(k = 0; k < NUM_ITERS; k += ARRAY_SIZE) {
		for(i = 0; i < ARRAY_SIZE; i++) {
			/* Calculate a[i+1] */
			mpf_ui_div(a[i+1], 1, aSqrt[i]);
			mpf_add(a[i+1], a[i+1], aSqrt[i]);
			mpf_div_ui(a[i+1], a[i+1], 2);
			mpf_sqrt(aSqrt[i+1], a[i+1]);

			/* Calculate b[i+1] */
			mpf_add_ui(b[i+1], b[i], 1);
			mpf_mul(b[i+1], b[i+1], aSqrt[i]);
			mpf_add(bBottom, a[i], b[i]);
			mpf_div(b[i+1], b[i+1], bBottom);

			/* Calculate p[i+1] */
			mpf_add_ui(pnp, a[i+1], 1);
			mpf_mul(pnp, pnp, pn);
			mpf_add_ui(pBottom, b[i+1], 1);
			mpf_mul(pnp, pnp, b[i+1]);
			mpf_div(pn, pnp, pBottom);
		}

		/* Recycle array: Last elements become first and process is restarted */
		mpf_set(a[0], a[ARRAY_SIZE]);
		mpf_set(aSqrt[0], aSqrt[ARRAY_SIZE]);
		mpf_set(b[0], b[ARRAY_SIZE]);
	}

	/* Print final results */
	gmp_printf("Final a[n]: %.6Ff\n", a[0]);
	gmp_printf("Final b[n]: %.6Ff\n", b[0]);
	gmp_printf("Final p[n] (PI): %.6Ff\n", pn);

	/* Free everything */
	for(i = 0; i <= ARRAY_SIZE; i++) {
		mpf_clear(a[i]);
		mpf_clear(aSqrt[i]);
		mpf_clear(b[i]);
	}
	mpf_clear(pn);
	mpf_clear(pnp);
	mpf_clear(bBottom);
	mpf_clear(pBottom);

	return 0;
}